A novel multi-agent simulation based particle swarm optimization algorithm

Recently, there has been considerable research on combining multi-agent simulation and particle swarm optimization in practice. However, most existing studies are limited to specific engineering fields or problems without summarizing a general and universal combination framework. Moreover, particle swarm optimization can be less effective in complex problems due to its weakness in balancing exploration and exploitation. Yet, it is not common to combine multi-agent simulation with improved versions of the algorithm. Therefore, this paper proposes an improved particle swarm optimization algorithm, introducing a multi-level structure and a competition mechanism to enhance exploration while balancing exploitation. The performance of the algorithm is tested by a set of comparison experiments. The results have verified its capability of converging to high-quality solutions at a fast rate while holding the swarm diversity. Further, a problem-independent simulation-optimization approach is proposed, which integrates the improved algorithm into multi-agent systems, aiming to simulate realistic scenarios dynamically and solve related optimization problems simultaneously. The approach is implemented in a response planning system to find optimal arrangements for response operations after the Sanchi oil spill accident. Results of the case study suggest that compared with the commonly-used shortest distance selection method, the proposed approach significantly shortens the overall response time, improves response efficiency, and mitigates environmental pollution.


Introduction
The merging of optimization and simulation technologies has seen rapid growth recently. The optimization of simulation models refers to finding the best values of some decision variables for a system where the performance is evaluated based on the output of a simulation model of this system [1]. However, due to the complexity, analytical simulation methods may be unsuitable for large-scale systems in manufacturing, supply chain management, financial management, etc. Agent-based modeling (ABM) is a bottom-up modeling method with high autonomy and interactivity, which is particularly suitable for complex system research [2]. Given the advantages of ABM, researchers have been solving optimization problems by integrating multi-agent simulation (MAS) with a wide range of methods, including game theory, reinforcement learning, and swarm intelligence (SI). Applications can be found in various fields such as intelligent transport systems, efficient electric grids, and smart buildings. Among the methods mentioned above, SI has been proved to be robust and efficient for optimization problems. As a classic category of SI, particle swarm optimization (PSO) has the advantage of simple parameter setting, fast convergence rate, high stability, and strong scalability, making it an excellent option for complex systems. However, despite a high convergence rate, the standard formation of PSO (SPSO) [3] is likely to stagnate at a local optimum, which reduces the solution accuracy. Thus, it is necessary to modify the algorithm for the application in complex problems where stable optimal solutions are indispensable. Different approaches have been considered to improve SPSO, among which setting parameters is especially representative [4]. Changing coefficients is a typical way of parameter setting. The performance of SPSO is affected by the values of its coefficients, i.e., the acceleration coefficients and the inertia weight. Thus, several attempts were made to tune the values of these parameters, such as the random weight method (RPSO) [5], the constriction factor approach (CPSO) [6], and the fuzzy adaptive approach (FAIPSO) [7]. Apart from the coefficients, the population size of the swarm is also a vital parameter. The basic idea for adjusting population size is to concentrate the individuals in the most promising area during the exploitation phase of the algorithm. De Oca et al. [8] conducted this idea by introducing incremental social learning (IPSO). On the other hand, as SPSO came from a simplified simulation of bird clustering behaviors, some researchers tried to improve it by mimicking other bird activities simultaneously. Neshat et al. [9] modified SPSO based on the predation phenomenon in bird swarms (PPSO). Besides these approaches, Li et al. [10] decoupled exploration and exploitation to make the algorithm more effective in large-scale optimization.
In MAS, agents can be either homogeneous or heterogeneous. Homogeneous agents are usually treated as particles when PSO is applied [11]. Yang et al. [12] presented a multi-agentbased model which simulated human behaviors in a multi-exit evacuation environment, where individuals were homogeneous agents and PSO was introduced to simulate their movement. Ali et al. [13] proposed a collective motion and self-organization control of a swarm of 10 unmanned aerial vehicles, which were divided into two clusters of five agents each. PSO was adopted to provide the best agents of the cluster. While homogeneous systems often simulate certain collective activities of agents with the same properties, heterogeneous systems are more practical for realistic scenarios due to the diversity of agents. Kanaga and Valarmathi [14] proposed a multi-agent model using PSO to solve patient scheduling problems, using a PSO Agent to perform the algorithm. Yu et al. [15] designed a Web services selection model, using PSO to select the optimal Web services combination in E-business. The model was realized by MAS. Thiel et al. [16] simulated Hanoi dwellers' choices between traditional markets and a hypothetical new one, aiming to find the optimal location for the ready-to-build local supermarket by using PSO to maximize sales volume. Mahad et al. [17] proposed an intelligent MAS approach to optimize power supply in community-based multi-microgrids systems, where various layers of autonomous and intelligent agents took decisions based on the PSO method.
Although there are many studies on improving SPSO, there are not so many attempts to combine proposed algorithms with MAS and apply them to practical problems. In fact, utilizing improved algorithms in simulation-optimization systems is necessary because complex tasks require high-performance optimizers, and SPSO may not meet this requirement. Besides, existing studies on the combination of MAS and PSO may vary in engineering fields, yet they all designed the system based on specific problems of their own research field without summarizing a general and universal framework, which is field-independent and can be implemented in various settings. Therefore, this paper focuses on a more in-depth exploration of improving PSO and integrating it into MAS to solve practical problems. The main contributions of this paper are as follows: 1. Proposing an improved particle swarm optimization algorithm (CoPSO) according to the biological nature of birds. It has a multi-level structure and a competition mechanism, aiming to improve performance by enhancing exploration while balancing exploitation.
2. Proposing a general simulation-optimization approach (MAS-CoPSO) to integrate CoPSO into MAS. It is supposed to achieve a very close combination between simulation and optimization, and the application should not be limited to certain engineering fields and specific problems.
3. Conducting a series of experiments to substantiate the validity of CoPSO.
4. Implementing MAS-CoPSO in a case study of the Sanchi oil spill accident to solve the response planning problem.

Methods and materials
This section elaborates on the methodologies used in this study. Firstly, a brief introduction to MAS technology is given. Detailed descriptions of the proposed CoPSO algorithm and the MAS-CoPSO approach follow afterward.

Multi-agent simulation
The fundamental element in MAS is the agents, which can be software or physical entities. Although different agents exhibit distinct behaviors, they share some common properties. For instance, they all have a certain degree of autonomy that enables them to work even without human intervention, they can communicate and interact with each other, and they are capable of perceiving and reacting to the changes in the environment as well as determining the proper behaviors to achieve the final goal [18]. In a simulation system, agents can be either homogeneous or heterogeneous depending on the specific problem background. Heterogeneous agents can deal with diverse scenarios and complex tasks, making it more suitable for complicated realistic systems [11]. Based on the aim of proposing a general framework for the combination of MAS and CoPSO, this paper mainly focuses on heterogeneous systems due to their universality.

Improved particle swarm optimization algorithm
PSO mimics the behavior of bird swarms. The fundamental idea is to assume a potential solution to the optimization problem as a bird without quality and volume, i.e., a particle. Particles are described by their positions and velocities. Generally, the position of a particle represents a specific solution while the velocity represents the searching direction and scope, which leads the particle to fly through the solution space. Every particle will get a fitness value by calculating the objective function. It will also adjust the value according to the experience of itself and the neighbors. Supposing the search space is D-dimensional, the ith particle is represented as , X d min and X d max are the lower and upper bounds of the dth dimension, respectively. The velocity of the ith particle is represented , V min and V max are the minimum velocity and maximum velocity specified by the user, respectively. In SPSO [3], during each iteration, particle i updates its position and velocity as follows: where w is the inertia weight, r 1 and r 2 are random values between 0 and 1, c 1 and c 2 are acceleration constants which control how the particle moves, is the best previous position of the ith particle, and gbest = (gbest 1 , gbest 2 , . . ., gbest D ) is the best position among all the particles in the swarm. PSO encourages exploration during the early stage of the iterations while encouraging exploitation in the latter stage. However, Angeline [19] has pointed out that the balance between exploration and exploitation is subtle and difficult to control in PSO, leading to stagnation or premature convergence. Inspired by the features of bird flocks, this paper proposes CoPSO to improve the performance of PSO, introducing a multi-level structure and a competition mechanism.
Multi-level structure. According to the biological nature of birds, their ability to sing and create a complete song is learned from their parents. Experiments have shown that if a bird is reared in silence, it can only scream. In addition, birds' singing skills are not set in stone. Basically, as they grow up, their expertise in singing also gets improved [20]. Based on this phenomenon, CoPSO has a multi-level structure, dividing the particles in a swarm into different singing levels. It is stipulated that the closer a particle is to a possible global optimum, the higher its level. The level of the ith particle is denoted by l i , which will increase if the particle doesn't update its personal best position pbest i in μ continuous iterations. Since birds of different ages learn to sing at different rates, particles of different levels have different acceleration constants. Fig 1 illustrates how the particles' levels change during the iterations in two-dimensional unimodal problems. At the early stage of the algorithm, particles are like newborn birds with little knowledge about singing. They will explore the search space and try to find promising areas. Once found, they might temporarily stop updating their personal best positions, leading to an upgrade. At the latter stage of the algorithm, some particles can be very close to the global optimum. Consequently, they will find it difficult to locate better positions, and their levels are also relatively higher.

PLOS ONE
Competition mechanism. While a particle's high level is likely to come from finding the global optimum, there are exceptions. Note that falling into local optima can also stop particles from updating their personal best positions. Since particles' levels are linked to their ages (the higher level and the more skillful in singing, the older), the idea of the competition mechanism inside bird swarms is utilized to improve the algorithm. The survival resources of swarms and the lifespan of birds are both limited. Therefore, a bird swarm is constantly renewed by the death of old birds and the born of new ones. If a particle's level is too high, it becomes too old to possess the resources of the swarm. The population size is set as a constant, once a particle's level reaches the upper bound L max , it will be replaced by a randomly generated new one. The new particle has the lowest level since it's a newborn. Fig 2 depicts how the competition mechanism functions in two-dimensional multimodal problems. If particles prematurely converge to a local optimum, they will upgrade around it. It can be seen from Fig 2a that few particle gets close to the global optimum. However, once the particles in the local area become too high-leveled, they will be replaced by new birds. As these new particles are randomly generated, they have the potential to explore unknown search space and finally locate the global optimum, which is shown in Fig 2b. Since particles tend to move closer to the global optimum and oscillate around it as the algorithm runs, in practice, the value of μ will be increased as the level rises, which allows an elaborate control over the balance between exploration and exploitation. The pseudo-code of CoPSO is shown in Algorithm 1.

Algorithm 1: Pseudo-code of CoPSO
Input: The value of L max , the value of μ, and the maximum number of iterations I max Output: The final solution gbest and its fitness value f(gbest) 1 Initialize gbest; 2 for each particle i do 3 Generate x i and v i randomly; is better than f(gbest) then 6 gbest = x i ; 7 end 8 end 9 while I max is not met do 10 for each particle i do 11 Update v i according to Eq (1); 12 Update x i according to Eq (2); if pbest i hasn't changed in μ continuous iterations then 20 Regenerate x i and v i randomly; 24

Problem-independent simulation-optimization approach
The architecture of the proposed MAS-CoPSO approach is shown in Fig 3, where the MAS module and the CoPSO module are the main components. For the optimization problem brought up by the user (e.g., decision making, scheduling, route planning, etc.), a suitable simulation time step should be chosen first according to specific backgrounds. In every time step, MAS is supposed to simulate the problem scene, during which the agents' behaviors might change the values of related variables and parameters. Simulation statistics are then put into

PLOS ONE
the optimization module, whose objective can be maximizing the efficiency or minimizing the cost of certain activities. CoPSO algorithm is supposed to optimize the objective function while satisfying all the constraints. The optimization result will then play as the input of MAS, instructing the simulation process in the next time step. This iteration will go on until meeting the stop criteria set by the user. The outcome of MAS-CoPSO is the final solution to the problem. For complex systems, the solution may not be strictly optimal (actually it is intractable to find the strict global optimum for these problems). However, MAS-CoPSO can compromise between the solution quality and computational cost, making it possible to get a satisfying result in a reasonable period of time.
Successful simulation of the problem scene relies on a proper setup of the environment (which reflects the background of the problem) and an appropriate design of all the agents. In MAS, agents have autonomous behaviors and complex interactions, and the relationship between them can be various (e.g., cooperation, competition, control, etc.). Value changes of related variables and parameters link simulation and optimization together. For CoPSO, simulation results will influence the calculation of the objective function's fitness value and the expression of constraints. Optimization results of CoPSO will function as the instructions for simulation in the next time step, i.e., how to adjust agents' behaviors and tune the values of variables, parameters, and environmental settings. These two modules compose a cohesive framework through consistent information exchange and close integration.

Experiments and analysis of the optimization algorithm
To substantiate the validity of the proposed CoPSO algorithm, this study conducts numerical comparison experiments and computational complexity analysis. Exhaustive descriptions and discussions are given in this section.

Comparison experiments
The performance of CoPSO is compared with a set of classic algorithms: • SPSO [3]: the standard formation of PSO, introduced in detail in the previous section.
• RPSO [5]: PSO with a random inertia weight factor designed for tracking dynamic systems.
• CPSO [6]: PSO with a new constriction factor related to the original acceleration coefficients.
• IPSO [8]: PSO with an increasing population size, i.e., whenever the algorithm can not find a satisfactory solution, add a new particle to the population.
Benchmarks and parameter setting. Five well-known benchmark functions commonly used in literature [21] are applied to evaluate the performance of CoPSO, both in terms of solution quality and convergence rate. They are non-linear minimization problems that present different difficulties to the optimizers. Table 1 lists the functions, the problem dimension n, the global minimum fitness value f min , and the search space ranges (which are also the initial ranges in this research).
The main parameters are set based on the recommendation in [22]. In CPSO, the maximal velocity V max and the minimal velocity V min are 100,000 and -100,000, respectively (since it is believed that V max and V min are not even needed in CPSO [6]). For other algorithms, V max is 10% of the search space's upper bound while V min is 10% of the lower one. The acceleration constants c 1 and c 2 are both 2.05 for CPSO. For SPSO, RPSO, and IPSO, the figure is 2.0. For CoPSO, L max is 3, the initial value of μ is 5 (added by 1 for each level up), and the gap between the acceleration constants of each level is 0.2 (c 1 and c 2 are the same with a maximum value of 2.0). A fixed inertia weight value of 0.6 is adopted for SPSO, IPSO, and CoPSO. The maximum iteration number is 1000 for all algorithms, meaning that the maximum particle number is 1000 in IPSO. In other cases, the population size is 80. A total of 50 repeating runs are conducted for each experiment. In addition, as recommended by Yang et al. [23], Wilcoxon rank sum tests are conducted for the statistical analysis between CoPSO and the other algorithms with a significance level of 0.05.
Results and discussion. Table 2 summarizes all the experimental results of 50 independent runs. Wilcoxon rank sum tests further verify the significance of these numerical results. The numbers in bold represent the comparatively best values. In general, CoPSO greatly outperforms SPSO, RPSO, CPSO, and IPSO on most metrics. For the exceptions, it still shows equivalent performances to the comparison algorithms. These results indicate that the multilevel structure and the competition mechanism can improve the solution quality of CoPSO.
It is believed that CoPSO can preserve excellent exploration and exploitation abilities while keeping a good balance between these two. In other words, CoPSO is able to jump out of local optima and fully exploit the promising areas at a fast rate. To substantiate this hypothesis, swarm diversity is used to measure an algorithm's exploration ability. Supposing the size of swarm S is N, the problem dimension is D, and x i represents the position of particle i. As recommended by Yang et al. [23], the diversity of swarm S is computed as follows: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where D(S) is the swarm diversity of S and � x is the average position of the swarm. Further, the converging speed can reflect whether an algorithm can compromise between exploration and exploitation properly. During the iterations, average swarm diversities and average global best fitness values (of 50 runs, in logarithmic form) for each algorithm on five benchmarks are recorded in Fig 4. Generally, CoPSO converges faster with better solutions than all other algorithms in most cases. As for swarm diversity, CoPSO holds a similar tendency on all benchmarks, i.e., decreasing rapidly in the early stage of the evolution and staying at a relatively high value (or even the highest) in the latter stage. This can be explained by the structure of CoPSO. Particles of different levels are in different regions of the search space. Those of lower levels are more likely to Schwefel's function The problem dimension n, the global optimum f min , and the search range are listed in the table. https://doi.org/10.1371/journal.pone.0275849.t001

PLOS ONE
be far away from the global best position while those of higher levels can be very close to it. CoPSO treats different levels differently by enhancing exploration for lower levels and promoting exploitation for higher ones. Thus, the algorithm has more potential to locate promising areas in the early stage and then converges to them very quickly, causing the swarm diversity to decline. However, particles may not upgrade because they find a global optimum, but because they are trapped in local optima. In this case, the competition mechanism will force the swarm to escape from local areas, which increases the diversity simultaneously. Despite the overall similarity, the behaviors of CoPSO differ slightly on different kinds of test functions. f 1 is a simple unimodal function with only one global minimum, making it possible to achieve fast convergence. Although all algorithms converge exponentially to the optima, CoPSO greatly surpasses the others due to better exploitation of the promising areas at the beginning and high-intensity exploration in the latter stage. For f 2 , CoPSO has a

PLOS ONE
competitive performance compared to RPSO and even shows better potential in the end as a result of higher swarm diversity. f 3 is a noisy quartic function. f 4 is a classic optimization problem with a global minimum inside a long, narrow, parabolic-shaped valley. f 5 is a multimodal function where the number of local optima increases exponentially as the problem dimension increases. For all these complicated benchmarks, CoPSO always keeps an excellent balance between exploration and exploitation, which helps it to achieve the fastest converging rate or the shortest stagnation at local optima. The other algorithms lack the adjusting ability of CoPSO. For instance, IPSO puts too much emphasis on exploration (holding the highest swarm diversity on most benchmarks) while RPSO and CPSO are too focused on exploitation, which degenerates their overall performances.

Computational complexity analysis
According to Yang et al. [23], given a fixed number of fitness evaluations, the computational complexity of an evolutionary algorithm is generally calculated by analyzing the extra cost in each generation without considering the cost of function evaluations, which is problem-dependent. As CoPSO inherits the simple structure of SPSO, its computational complexity will be analyzed by comparison with SPSO.

PLOS ONE
possible in certain generations. Besides, the value increase of μ also controls the time cost of this part. As for the space complexity, since CoPSO needs to store the information of l i and the time particles stay at the same pbest i , it requires O(N) extra memory than SPSO (which takes O(N × D) space to save particles' positions and velocities). To conclude, with proper parameter setting, the computational complexity of CoPSO can be controlled within an acceptable range. Table 3 lists the average computing time for a single run (in 50 repeating experiments) with the same parameters as the numerical experiments. SPSO, RPSO, and CPSO take nearly the same time on different functions due to their similar algorithmic structures. However, IPSO has a growing population, dramatically increasing fitness evaluations at the latter stage of the iterations. Consequently, the computational cost of IPSO is the highest. Based on the complexity analysis above, CoPSO takes acceptable extra time in each iteration compared to SPSO, verified by the experimental result that CoPSO costs slightly higher than SPSO on the benchmarks. In conclusion, CoPSO surpasses the other algorithms in solution quality while remaining computationally efficacious.

Case study of the simulation-optimization approach
To demonstrate that MAS-CoPSO can be used in practical applications, a case study of a real marine oil spill accident is carried out. This section presents the accident background, the implementation of MAS-CoPSO in a response planning system, and the mathematical formulation of the optimization problem. Further, an in-depth analysis of the system's performance is given.

Background
Frequent marine oil spills have posed severe tests to the environment. Reducing the loss caused by oil spills has become an important issue. The capability of the MAS-CoPSO approach is tested through a case study of the Sanchi oil spill accident, which is caused by a collision between Panamanian oil tanker Sanchi and Chinese bulk carrier Changfeng Crystal. The accident happened on January 6, 2018, and the collision site is about 160 nautical miles east of the Yangtze River Estuary. Fig 5 illustrates the approximate location where the accident happened and the areas around it. Tons of condensate oil carried by Sanchi leaked into the ocean, causing significant economic loss and environmental pollution. Researchers have pointed out that oil tankers are transporting around 90% of all the oil around the world [24]. Furthermore, ship collision is the primary cause of many catastrophic marine oil spills, which suggests the representativeness of the case study.
Zhong and You [25] have concluded that after spill accidents, oil leaked into the ocean forms scattered slicks around the spill site and undergoes various physical and chemical changes due to environmental forces. The most dominant processes are spreading, evaporation, emulsification, and dispersion. These processes will dramatically change the volume, area, thickness, and viscosity of those slicks, which will affect the efficiency of oil spill response operations as a consequence.
There are some equipped response teams berthed at docks near the spill site. They will conduct a wide range of actions to clean up the ocean, among which booms, skimmers, chemical dispersants, and in situ burning are the most frequently used methods. Booms are generally the first equipment deployed after an oil spill and are often used to protect shorelines, divert oil to certain areas, or concentrate oil for further recovery and burning. Skimmers are adopted to recover oil or oil-water mixtures from the water surface. The characters of the slick will affect the working efficiency of skimmers. Chemical dispersants can reduce the oil-water interfacial tension, which should be initialized as soon as possible. In situ burning indicates controlled burning of the oil at or near the site. In this case, it is assumed that the spraying of dispersants and the concentration of leaked oil have been completed quickly after the accident. Besides, in situ burning will emit toxic substances and waste the leaked oil. Thus, the primary response operation is recovery.
This study focuses on supporting response planning after the accident. A simulation-optimization system is built based on the proposed MAS-CoPSO method. The overall objective is to clean up the polluted sea area in the most efficient way, i.e., as soon as possible. The system is supposed to realize interactive spill simulation and response optimization while considering the influences of the oceanic environment. It should also give the most reasonable allocations for available resources and the schedule of response teams' actions.

Implementation
The implementation of MAS-CoPSO in the system is described in Fig 6, which is an instantiation of Fig 3. The MAS module and the CoPSO module work collectively, composing a dynamic system for marine oil spill response planning. The CoPSO algorithm is encapsulated as a function to be called before each time step. Meanwhile, agents of the simulation module consistently update their behaviors and interact with each other, realizing the oil spill fate modeling and response operation modeling. The agents must follow specific rules in reality. For instance, the property changes of oil slicks are calculated by weathering models considering spreading, evaporation, emulsification, and dispersion. Skimmers and ships of response teams should obey the rules for oil recovery. The performance of a response team can be affected by the actions of other teams and the characteristics of oil slicks. For example, when a skimmer works on a slick to collect oil, the volume and area of the slick, the evaporated and dispersed oil rates, the viscosity and the water content of the oil, will all be affected. If more than one team are instructed to head for the same slick, they will collaborate.
To complete the global objective of minimizing response time while cleaning up the spill site to a large extent (the stop criteria, only a small portion of leaked oil remained in the ocean), CoPSO is utilized to maximize the decreased volume of oil in each simulation time step. Meanwhile, MAS should simulate equipment location, response process, oil spill fate, and transportation. The system can be further updated to suit different purposes and requirements.

Mathematical formulation of the optimization problem
Since the weathering process dramatically affects oil properties, it is necessary to simulate it as accurately as possible. Furthermore, the objective function and constraints of the optimization problem can be derived from the model.
Simulation of the weathering process. First of all, the initial area of a slick, A 0 (m 2 ), can be calculated by the equation shown as follows [26]: where ρ ω (g � cm −3 ) is the density of seawater, ρ 0 (g � cm −3 ) is the density of oil, ν ω (0.801 × 10 −6 m 2 � s −1 under 30˚C) is the kinematic viscosity of seawater, V 0 (m 3 ) is the inital volume of the slick, g (m � s −1 ) is the acceleration of gravity, k 2 and k 3 are constants with values of 1.21 and 1.53, respectively. The area changing rate of the slick due to spreading can be modeled as follows [27]: where A (m 2 ) is the current surface area of the slick, t (s) refers to the time since the accident happened, V (m 3 ) is the current volume of the slick, K 1 (s −1 ) is a dominant physicochemical parameter of the oil with a default value of 150.

PLOS ONE
Evaporation is assumed as the most influential environmental effect, causing a loss of up to 20-50% of the spill. The rate that oil evaporates from the sea surface is modeled by the following equation [28]: where F E (%) is the current volume fraction of the oil that has been evaporated with an initial value of F E(t=0) = 0, T K (K) is the oil temperature, A ev and B ev are empirical constants with fixed values of 6.3 and 10.3, respectively. K ev (m � s −1 ) is the mass transfer coefficient for evaporation which can be calculated as follows [29]: where v wind (m � s −1 ) is the wind speed. T O and T G are the initial boiling point and the gradient of the oil distillation curve, respectively. Their values can be calculated through functions of the oil API (American Petroleum Institute) degree as follows: T O ¼ 457:16 À 3:3447 � API; T G ¼ 1356:7 À 247:36 � ln API: Emulsification can also influence the slick characteristics dramatically, especially for the oil viscosity. The initial oil viscosity μ 0 (%) can be calculated using the following equation [29]: where AC (%) is the asphaltene content of the parent oil, which is a constant. As the viscosity changes over time, its changing rate is given by [30]: where μ (%) is the current viscosity, C 3 is a constant for the final water content, C 4 is an oildependent constant, Y W (%) is the fractional water content in the emulsion, it changes over time, which can be computed with the following equation [27]: where K em is an empirical constant between 1 × 10 −6 and 2 × 10 −6 , Y W has an initial value of Y W(t = 0) = 0. Finally, dispersion is also a vital factor in the loss of oil. The volume of oil naturally dispersed changes over time, which can be modeled by [27]: where V D (m 3 ) is the total dispersed oil volume until now with an initial value of V D(t=0) = 0 and z t (cP) is the oil-water interfacial tension. Objective function and constraints. In the response process, the optimization problem is formulated on the basis of each simulation time step. Suppose that the ith response team is recovering oil on the kth slick during time step T step , define ORR i (m 3 � hr −1 ) as the oil recovery rate of team i (i.e., the amount of oil that the team can recover per hour) and ST k (m) as the current oil thickness of slick k. ORR i is determined by the properties of the team's skimmer and ST k : where α and β are empirical coefficients obtained from experimental tests, reflecting the features of the skimmer. ST k can be calculated by: where V k (m 3 ) and A k (m 2 ) are the current oil volume and the area of the kth slick, respectively. A k can be calculated according to Eqs (3) and (4), V k can be computed as follows: where V k0 (m 3 ) is the remaining oil volume of this slick at the end of the last time step (i.e., the initial oil volume of the current time step), EV k (m 3 ) is the volume loss caused by evaporation, DV k (m 3 ) is the volume loss caused by dispersion, and SV i (m 3 ) is the volume recovered by team i. Notice that SV i may not exist, which indicates that no response team is working on slick k right now. The system can get the value of EV k , DV k , and SV i according to Eqs (5), (8) and (9), respectively. During the simulation, for simplification, a short T step should be chosen and the system will use the value of V k0 for the calculations regarding V k . The aim of the system is to clean up the polluted sea area as soon as possible. Therefore, based on the proposed method, given a fixed simulation time step, the oil volume loss will be maximized during each step, which is composed of the evaporation loss, the dispersion loss and the recovered volume. Suppose there are m oil slicks and n response teams during T step , V step represents the total oil volume loss in this period. F Ek and μ k refer to the evaporated oil fraction and the oil viscosity of slick k, respectively. The objective function and constraints are given by: s.t.

Performance of the system
To examine the efficiency of MAS-CoPSO, its performance is compared with the shortest distance selection approach (SDS). Ye et al. [31] have suggested that SDS is a simple strategy commonly used in marine oil spill emergency response, which indicates a process that allows a response team to choose the nearest oil slick as the target for oil recovery. After meeting the cleaning requirements, they will choose the second nearest oil slicks to continue. Therefore, when conducting SDS, the distances between ships and slicks are the judgment criteria for response planning. When initializing the simulation process, each response team (corresponding to a certain response agent) is assigned with the closest oil slick (corresponding to a certain spill agent) to clean. Before each time step, instead of calling the CoPSO module, SDS will check whether the slick is cleaned up, if the oil volume is below the set threshold, the ship will head to the closest one in the remaining slicks, otherwise, it will stay at the current location until finishing the cleaning task. Compared to MAS-CoPSO, SDS is more straightforward, which explains why it is commonly used in real accidents. However, SDS makes plans without considering the oceanic forces and the interactions between agents, which according to previous discussions, can dramatically influence the overall efficiency. Thus, comparative experiments are conducted for further analysis.
Simulation settings. Chen et al. [32] introduced the process and consequences of the accident in detail. Based on their description, it is assumed that the spilled oil was split into ten slicks within the East Sea with a total volume of 113,000 tons. Table 4 lists the oil volumes and locations (described by the distance and direction from the approximate collision site) of these slicks. The ship-mounted skimmers belonging to five different response teams are the only available nearby cleanup means to be applied. Assume that all the teams have enough storage space for recovered oil. As the allocation process needs specific transportation time, the speed of ships is set at 30 km � h −1 . Skimmers of different teams have different recovering efficiencies, which are affected by α and β in Eq (9). Table 5 lists the values of the empirical coefficients for different teams. Fig 7 illustrates the simulation process. The red cross represents the collision site, the black spills represent the oil slicks (the size reflects the oil volume of the corresponding slick), and the ships represent the response teams carrying the recovery task. The position and size of a slick will change over time. A ship can either be heading to a slick or working on it.
As mentioned above, spreading, evaporation, emulsification, and dispersion are the main oceanic effects. Table 6 lists the inputs for the modeling of these processes, including

PLOS ONE
environmental parameters and the characteristics of the leaked oil (condensate oil). As [26] suggested, the drifting speed of oil slicks is about 2.5-4.5% of the wind speed, so it is set at 0.03 m � s −1 .
The model is built in AnyLogic1. 50 repeating runs are carried out with 50 iterations per time step (T step is one hour). Once the volume of the remaining oil is reduced to 10% of the original figure, the simulation will stop.
Results and discussion. Table 7 presents the simulation results of MAS-CoPSO and SDS with the same stop criteria and time step. The operation time for achieving an oil recovery rate

PLOS ONE
of 90% is 129.82 hours based on the optimal outcomes given by the MAS-CoPSO system, which is only 74.2% of the figure for SDS. Shorter response time means higher operation efficiency, less volume of wasted oil, and slighter damage to the environment. The results also indicate that total recovered oil holds a proportion of 65.54% in MAS-CoPSO, which is approximately 1.2 times higher than the number for SDS. Fig 8 illustrates the variations of remaining oil volume for each slick in MAS-CoPSO and SDS scenarios during the entire simulation process. The optimization outcomes of MAS-CoPSO are highly related to the thickness of the oil slicks, and the system intends to keep a balanced volume level for each slick by optimizing the time for allocation and the response efficiency (see Fig 8a). However, oil thickness does not affect the decision-making of the SDS model, so it will not try to reduce the volume steadily. As a result of different strategies, the curves of MAS-CoPSO are much smoother than those of SDS. In addition, the SDS method might ignore slicks too far from the response teams (e.g., Slick 1), leading to long-time exposure of the oil, emitting more harmful substances (see Fig 8b). Fig 9 depicts the variations of accumulated recovered volume for each team in two approaches. MAS-CoPSO tries to find optimal schedules for all teams in every time step, so the curves climbs stably (see Fig 9a). Meanwhile, SDS always cleans up one slick before moving to another, leading to discontinuous curves, which also hinders the increase of the overall efficiency (see Fig 9b). The operation time (hr), the oil volume recovered by each team (ton), the proportion of each oil loss type (%), and the proportion of remaining oil after the response process (%) are listed in the table. https://doi.org/10.1371/journal.pone.0275849.t007

PLOS ONE
The above results show that the response planning system based on MAS-CoPSO outperforms the commonly used SDS strategy in terms of time consumption, oil recovery rate, and environmental impact. Therefore, oil spill response teams can schedule their operations according to the simulation-optimization outcomes of this system rather than relying on the inefficient SDS approach. Moreover, complex problems and high-intensity interactions can enhance the advantages of MAS-CoPSO, indicating the system's ability to function in scenarios with higher oil volumes and more response teams. Even though the case study is simplified and focuses on the oil recovery process, the system has the potential to comprehensively support multiple cleanup techniques concerning booms, chemical dispersants, and in situ burning. Besides, modeling of the weathering process can be easily enriched by considering more complicated oceanic forces such as dissolution, photo-oxidation, sedimentation, and biodegradation [25]. Hydrodynamic simulation of oil spill trajectories can also be considered. In addition, the application range of the MAS-CoPSO-based system can be further enlarged by handling uncertainties and risk assessments. With these extensions, oil spill response teams can get a practical simulation-optimization tool to support their planning process.

Conclusion
An improved particle swarm optimization algorithm called CoPSO is proposed in this paper. Inspired by the features of bird flocks, a multi-level design and a competition mechanism are introduced into CoPSO to balance the exploration and exploitation more effectively, increase solution accuracy, achieve fast convergence, and avoid stagnation at local optima. The performance of CoPSO is tested on a series of well-known benchmark functions. Comparisons are made between CoPSO, SPSO, and three classic variants of it. Experimental results have verified its capability.
This paper also proposes a problem-independent simulation-optimization approach called MAS-CoPSO to combine CoPSO with MAS. MAS-CoPSO can achieve a dynamic simulation of various complex systems and give satisfying solutions to related optimization problems. In MAS-CoPSO, the relationship between the simulation module and the optimization module is highly cohesive. The effectiveness of MAS-CoPSO is substantiated through a response planning system for the Sanchi oil spill accident. The system simulates the accident scene with consideration of the oil slicks' physicochemical evolution. Case study results show that the system can optimize response device allocation, operation scheduling, and time consumption, causing slighter damage to the environment.
MAS-CoPSO can also be applied in other scenarios to solve different optimization problems depending on the users' requirements. For future studies, we will consider the combination with other methods, such as fuzzy control and neural network.